Impacts of ocean warming on echinoderms: A meta‐analysis

Abstract Rising ocean temperatures are threatening marine species and populations worldwide, and ectothermic taxa are particularly vulnerable. Echinoderms are an ecologically important phylum of marine ectotherms and shifts in their population dynamics can have profound impacts on the marine environment. The effects of warming on echinoderms are highly variable across controlled laboratory‐based studies. Accordingly, synthesis of these studies will facilitate the better understanding of broad patterns in responses of echinoderms to ocean warming. Herein, a meta‐analysis incorporating the results of 85 studies (710 individual responses) is presented, exploring the effects of warming on various performance predictors. The mean responses of echinoderms to all magnitudes of warming were compared across multiple biological responses, ontogenetic life stages, taxonomic classes, and regions, facilitated by multivariate linear mixed effects models. Further models were conducted, which only incorporated responses to warming greater than the projected end‐of‐century mean annual temperatures at the collection sites. This meta‐analysis provides evidence that ocean warming will generally accelerate metabolic rate (+32%) and reduce survival (−35%) in echinoderms, and echinoderms from subtropical (−9%) and tropical (−8%) regions will be the most vulnerable. The relatively high vulnerability of echinoderm larvae to warming (−20%) indicates that this life stage may be a significant developmental bottleneck in the near‐future, likely reducing successful recruitment into populations. Furthermore, asteroids appear to be the class of echinoderms that are most negatively affected by elevated temperature (−30%). When considering only responses to magnitudes of warming representative of end‐of‐century climate change projections, the negative impacts on asteroids, tropical species and juveniles were exacerbated (−51%, −34% and −40% respectively). The results of these analyses enable better predictions of how keystone and invasive echinoderm species may perform in a warmer ocean, and the possible consequences for populations, communities and ecosystems.


| INTRODUC TI ON
Ocean warming, a repercussion of anthropogenic climate change, is one of the principal emerging threats facing marine ecosystems (IPCC, 2019). An increased prevalence and intensity of marine heatwaves (i.e. warming above the 90th percentile for at least five consecutive days; Hobday et al., 2016) in recent years has already resulted in the degradation of some of the world's most vulnerable marine ecosystems, such as coral reefs, sea grass beds and kelp forests (Filbee-Dexter et al., 2020;Mellin et al., 2019;Strydom et al., 2020). Within these ecosystems, ectotherms are particularly thermosensitive because their body temperature, and therefore their physiology and behaviour, are naturally linked to their thermal environment (Lagerspetz & Vainio, 2006). Accordingly, ocean warming has been implicated in population declines, phenological shifts and range shifts in an array of marine ectotherms (Poloczanska et al., 2016). Considering the ecological and economic importance of many marine ectotherms Ling et al., 2009;Munday et al., 2013), a knowledge of how elevated temperatures will impact these organisms is essential.
Due to the widespread threat of ocean warming, it is important to understand how common biological and ecological aspects shape the responses of marine ectotherms to a warming environment. For example, certain biological responses (e.g. metabolic activity and survival rates) and life stages may have differing sensitivities to elevated temperatures (Dahlke et al., 2020;Pereira Santos et al., 2021;Sampaio et al., 2021). Biological responses associated with physiology are likely to be particularly sensitive to environmental warming in ectotherms, with their magnitude and direction dependent upon where the temperatures experienced fall within the organisms' thermal window (Pörtner et al., 2017;Sampaio et al., 2021). Physiological performance is often assessed nonlethally via oxygen consumption rate measurements, which serves as a proxy for aerobic metabolism, and evidence suggests that this measure can operate as a good proxy for other aspects of performance (Pörtner et al., 2017;Svendsen et al., 2016). In general, aerobic metabolism is only able to increase to the organisms' thermal limit, after which, constraints on oxygen delivery will lead to metabolic depression and declines in energy production (Lang et al., 2021;Pörtner et al., 2017). The additional aerobic metabolic costs of warming may impact growth, development, activity, feeding and reproduction (Pörtner et al., 2017;Schulte, 2015).
Reproduction is typically sensitive to elevated temperatures, and resulting embryos and larvae are often highly vulnerable to warming Collin et al., 2021;Dahlke et al., 2020), and other concurrent stressors, such as ocean acidification Kroeker et al., 2010;Pandori & Sorte, 2018;Przeslawski et al., 2015). On the contrary, juveniles and adults are usually more thermotolerant, testament to their well-developed physiological capacity and greater ability to modify their behaviour to withstand warmer conditions Dahlke et al., 2020;Lang et al., 2021;Nguyen et al., 2011;Sampaio et al., 2021). Behavioural plasticity and thermoregulation could allow marine ectotherms to buffer their experience of stressful thermal conditions and potentially reduce the costs of increased metabolic demands in warmer oceans (Wong & Candolin, 2015).
Responses of marine ectotherms to warming may be in part, a reflection of their evolutionary and thermal history. Responses to elevated temperatures in marine ectotherms have been demonstrated to vary taxonomically (e.g. Nguyen et al., 2011;Sampaio et al., 2021). These patterns are due to taxa having differing capacities to adjust their physiology and behaviour in response to warming, which may be associated with their phylogeny and thermal experience throughout evolutionary time (Bennett et al., 2021;Nguyen et al., 2011;Sunday et al., 2011). Tolerance to elevated temperatures is often attributed to the amount of thermal variation experienced (Hughes et al., 2018;Sunday et al., 2011). Seasonal thermal variability is low at tropical latitudes (Hughes et al., 2018;Madeira et al., 2017;Sunday et al., 2011). Consequently, marine ectotherms living in these relatively stable thermal environments often have a narrow thermal breadth, and live closer to their upper thermal limits compared with those inhabiting more variable thermal environments (Poloczanska et al., 2016;Pinsky et al., 2019;Schulte et al., 2011;Sherman, 2015;Sunday et al., 2011;Woolsey et al., 2015). Marine ectotherms exposed to greater seasonal thermal variability at temperate latitudes for instance often have wider thermal windows and have evolved greater capacity to withstand periods of thermal extremes via enhanced physiological and behavioural thermoregulation (Madeira et al., 2017;Sokolova & Pörtner, 2003;Sunday et al., 2011;Woolsey et al., 2015). This geographic pattern in vulnerability to warming may not always hold true, however, because vulnerability may also be dependent on other spatial factors, including the breadth of the geographic range, the part of the range in which they occur, their depth and habitat (Drake et al., 2017;O'Connor et al., 2012;Pey et al., 2011;Poloczanska et al., 2016;Sasaki et al., 2022;Zettlemoyer & Peterson, 2021).

T A X O N O M Y C L A S S I F I C A T I O N
Behavioural ecology, Biogeography, Ecophysiology, Evolutionary ecology, Global change ecology, Global ecology, Life history ecology, Population ecology have resulted in significant coral loss on Australia's Great Barrier Reef (Mellin et al., 2019), while shifts in the distribution of large populations of the echinoid Centrostephanus rodgersii have contributed to the widespread and accelerated loss of highly diverse kelp forest ecosystems in Tasmania (Ling, 2013). Ocean warming is considered to be the primary cause of the distributional shift of C. rodgersii (Ling et al., 2009). Ocean warming has also been implicated in increasing the incidence and severity of disease outbreaks in a variety of echinoderm species (Clemente et al., 2014;Dungan et al., 1982;Lester et al., 2007;Menge et al., 2016). For instance, the asteroid Pisaster ochraceus, the first species to be coined a 'keystone' species (Paine, 1966;Wagner, 2010), helps maintain invertebrate diversity by providing top-down control on mussel populations along the Pacific coast of North America, but has recently experienced significant population declines due to Sea Star Wasting Disease, which may be linked to ocean warming (Menge et al., 2016). Considering the numerous ecosystem impacts resulting from a 'boom' or 'bust' in echinoderm populations (Menge et al., 2016;Uthicke et al., 2009), there is a strong impetus to understand how echinoderm species may be impacted by a changing climate.
Despite an increasing number of controlled laboratory-based studies focussing on understanding the effects of elevated temperatures on echinoderm species, there is little synthesis of the overarching impacts of warming, or explanations for the variation in vulnerability to elevated temperature observed within this phylum. Meta-analyses have provided a useful statistical tool for quantitatively assessing broad patterns in the responses of marine ectotherms such as fish, molluscs and crustaceans to environmental change, including ocean warming, ocean acidification, hypoxia and reduced salinity (e.g. Hu et al., 2022;Kroeker et al., 2010;Pandori & Sorte, 2018;Pereira Santos et al., 2021;Przeslawski et al., 2015;Sampaio et al., 2021). To the authors' knowledge, there are no comparable meta-analyses to date that have focussed on understanding whether echinoderms adhere to general biological, ecological or evolutionary expectations regarding their vulnerability to ocean warming Dahlke et al., 2020;Poloczanska et al., 2016;Sampaio et al., 2021;Sunday et al., 2011).
Herein, a meta-analysis using a comprehensive empirical dataset (see File S1 for the full dataset) is presented, consisting of information from 85 studies that explore the effect of warming on echinoderms. Specifically, the direction and magnitude of the effect of warming (relative to control temperatures) on performance was investigated, and how this varies depending on the (1) biological response tested, (2) ontogenetic life stage, (3) taxonomic class, (4) region and (5) experimental design. To further explore the potential effect of future ocean warming, a subset of the data was examined in which responses were only included if the experimental temperature was above the predicted end-of-century mean annual temperature (MAT; IPCC, 2019) at the location of collection.
The outcomes of this meta-analysis will aid in assessing the vulnerability of key echinoderm species to current and future ocean warming, and the broader impacts on the marine environment.

| Literature search
Relevant publications on the effect of warming on echinoderms were initially identified using systematic searches within ISI Web of Science (studies published before June 2021). The keywords used were: 'temp*', 'thermal', 'thermotolerance', 'warm*' and 'climate change' (an asterisk represents the wildcard operator in Web of Science, used to find different endings of keywords). These words were combined (using the 'AND' boolean) with both the common and/or Latin names of the five classes of echinoderms, that is, Asteroidea (starfish/sea stars), Echinoidea (sea urchins and sand dollars), Ophiuroidea (brittle stars), Holothuroidea (sea cucumbers) and Crinoidea (feather stars and sea lilies). The latter class was not included in this review due to a lack of relevant studies. The field tag used was 'TS = Topic'. To avoid any bias associated with conducting literature searches in a single database (Martín-Martín et al., 2018), systematic searches using ISI Web of Science were supplemented using Google Scholar and by exploring studies referenced in relevant publications. There was no limit on the year of publication for included studies, but field-based studies were excluded, due to limited ability to control for other factors that may influence responses.

| Data collection and selection
The papers that met these initial criteria (n = 143) were further screened for eligibility (see Appendix 1 for a flow diagram of the selection process). Studies were only incorporated if the responses fell under one of the following key biological responses: developmental success, feeding and nutrition, growth, metabolic rate (routine or resting), movement, reproductive success and survival (n = 3 studies that did not fall under the aforementioned biological responses were therefore excluded). Acute 'ramping' (i.e. warming rate > 2°C h −1 ) studies were not included due to the associated cumulative stress likely causing confounding deleterious effects on organism responses to warming (n = 5 studies excluded). Studies and individual responses were excluded if the control temperature was not specified or clearly stated (n = 21 studies excluded), or if the mean response for a control and at least one experimental temperature were not reported (n = 4 studies excluded). Finally, studies were excluded if the standard deviation or sufficient information for it to be calculated/estimated (i.e. standard error, upper and lower 95% confidence limits, interquartile range, minimum and maximum values) and/or the sample sizes were not provided (n = 22 studies excluded). If multiple synonymous responses were included in the publication (e.g. gonad wet weight and gonad dry weight), only one was included in the analyses. The responses at all elevated temperatures were included regardless of how realistic these temperatures were in the context of near-future climate change. If a response was measured at multiple time points, in general the response at the first and last time point were used in the analyses, to avoid pseudoreplication, and taxonomic biases. If the first time point was immediately after environmental conditions were changed, the second and last time point were used instead (see Appendix 2 for a reference list of the data sources for the analyses).
For each response, the mean, standard deviation and number of replicates at both the control and experimental (elevated) temperatures were recorded. The control temperature was the temperature at which the individuals were habituated to prior to the experiments, or if there was no laboratory habituation period, then the temperature at the collection location was used. When raw or summary data tables were not published, but were presented in a graphical format, data were mined from the primary literature using the program WebPlotDigitizer (Automeris.io, 2021). The biological response being tested, ontogenetic life stage, the taxonomic class of the test subjects and the region were also recorded (see Appendix 3 for a table of these predictors and the groups within them). An attempt was made to gather data on the collection depth, and the latitudinal range of each species (to consider how range size, and the position of the individuals within the range, may impact vulnerability to warming), however, insufficient data were available (Living Australia, 2021). Experimental variables were recorded or calculated for each study, namely the habituation time at ambient temperature prior to experiments (days), the warming rate from the control to the experimental temperature (°C h −1 ), the exposure time at experimental temperatures (days) and the natural logarithm of the ratio of the experimental and control temperatures (LnSR; Pereira Sampaio et al., 2021).
The resulting values for LnSR were dependent not only on the difference between the temperature values, but also depended upon how hot or cold the control values were themselves; thus, the temperatures were adjusted accordingly. As in Pereira Santos et al. (2021) and Sampaio et al. (2021), all control temperatures were set to 2°C, and new experimental temperatures were calculated as the original experimental temperature minus the original control temperature, plus the new baseline control temperature.
The mean annual temperature (MAT) for each study location was also established by extracting COBE long-term mean sea surface temperature estimates from the National Oceanic and Atmospheric Administration (NOAA Physical Sciences Laboratory, 2021) using the coordinates of the collection locations (or approximate coordinates if not available). For each study, the temperature data from the 10 years prior to the year of the experiment were used to establish the MAT.

| Effect size and variance calculation
Effect sizes, that is the ln-transformed response ratios (LnRR), were calculated using the equation outlined in Hedges et al. (1999): where X E and X C are the mean responses in the experimental and control treatments, respectively. The natural logarithm linearises these values, causing the values of both the experimental and control means to be weighted equally, and removes some of the skewness of the sampling distribution (Hedges et al., 1999). A positive LnRR indicates that the response is positively impacted by elevated temperature, while a negative LnRR means that the response is negatively impacted.
Furthermore, if the LnRR is zero, then there is no effect of elevated temperature on the response. If for a given response, a higher mean value indicated a more detrimental effect (i.e. per cent abnormality, per cent mortality, development time and righting time), the sign of the LnRR was reversed for a more intuitive visualisation. The variance (inverse-variance weights) of each response was calculated using the equation of Hedges et al. (1999): where, S is the standard deviation and n is the sample size. The variance enables the precision of the estimate to be established. Observations with a greater sample size and lower standard deviation are more heavily weighted, as they are a more precise estimate of the effect size (Kroeker et al., 2010).

| Data analyses
Multivariate multilevel linear mixed effects models were conducted in R v.4.1.2 (R Core Team, 2021) and used to establish the mean effects of various predictor variables (see Appendix 3). Separate models were fitted for each predictor (i.e. biological response, life stage, taxonomic class, region, habituation time at the control temperature, warming rate, and exposure time at the experimental temperature), as there were insufficient data available to consider multiple predictors simultaneously. Responses for metabolic rate were only included in the analysis for the 'biological response' predictor, due to the uncertainty regarding whether positive or negative changes in the metabolic rate translate to beneficial or detrimental effects on organismal fitness (Hu et al., 2022;Pörtner et al., 2017;Schulte, 2015). All models were carried out using the rma.mv function from the 'metafor' package in R (Viechtbauer, 2010). These models account for the hierarchical structure of meta-analytic data and considers the variation both within and between studies, while also accounting for any nonindependence of effect sizes (Cheung, 2019;Jackson et al., 2011;Konstantopoulos, 2011).
The formula for each of the models, adapted from Pereira Santos et al. (2021) and Sampaio et al. (2021), was as follows: The categorical predictors (i.e. biological response, life stage, taxonomic class and region), were interacted with LnSR ('LnSR: Predictor -1'), which allowed for varying slopes of LnSR on LnRR for each predictor group, with a fixed intercept of 0 as in previous studies . LnSR was not required in the experimental variable models, as all predictors were continuous. The variance (v) was included in all models to weight responses based on their precision, as well as the random effects structure. 'Response number' was the number allocated to all responses at the different experimental temperatures for a given measure (e.g. per cent fertilisation success, gonad index, wet weight and locomotion speed) within a single study (see Appendix 4 for a full list of measures included within each 'biological response' group). The random effect was nested within the study number (~1|Study number/Response number). In most cases, there was only a single study per publication; however, if the paper included independent experiments on multiple echinoderm species, these were included as multiple studies per publication. Restricted maximum likelihood was used to fit all models, and t-statistic methods were implemented, which are more conservative than default z-statistic methods (Knapp & Hartung, 2003 and Strongylocentrotus droebachiensis (temperate and polar echinoid; n = 53 data points) were also tested in the taxonomic class and region models and were found not to substantially drive the results (see for methods and results of the sensitivity, publication bias and linearity tests).

| End-of-century scenario
All the predictor models above were re-run using a subset of the data where responses were only included if the experimental temperatures were above the MAT predicted for the end of the century (+2.58°C, RCP8.5 global mean sea surface temperature projection 2081-2100; IPCC, 2019) at the collection sites. This allowed the assessment of whether there is evidence of an increased negative impact on echinoderms when experiments use magnitudes of warming tantamount to end-of-century predictions, and whether testing temperatures that are already frequently experienced by the organism may underestimate the potential effects of future warming (see File S2 for code for all models, including sensitivity analyses and publication bias detection).

| RE SULTS
A total of 85 studies met the selection criteria for the 'biological   Figure 2c,g). Asteroidea were the most negatively affected, with an average 30% decline in performance (p < .001).
Ophiuroidea appeared to respond positively to thermal challenge, with a 15% increase in performance, but this was not significant (p = .244). Considering the relatively even spread of data points around zero for this class, the high estimate may be a consequence of low resolution due to the small sample size (n = 20 responses). For Echinoidea, the mean LnRR was close to zero (p = .638). Removal of the three species that contributed the most data points changed the model output very little (see Appendix 6.1).

| End-of-century scenario
When only the responses that applied experimental temperatures relevant to end-of-century projections were tested, it was found that for some of the groups within predictors, the effect size and/ or significance had changed (see Appendix 7 for a figure illustrating these results). For instance, metabolic rate was no longer significantly positively impacted by warming, and gametes were no longer significantly negatively impacted by warming (p > .050), yet the negative impact of warming on juveniles increased by 31% (p < .001).
Under this end-of-century scenario, Asteroidea and Holothuroidea were also more negatively affected, with performance declining . The mean effect sizes and the robust 95% confidence intervals are provided (left panels, a-d), as well as the raw data points (right panels, e-h). The number of responses in each group are provided in parentheses in the right panels. Significant mean effect sizes are indicated by an asterisk (*p < .05; **p < .01; ***p < .001) in the left panels. Mean effect sizes and 95% confidence intervals have been corrected for differences between the control and experimental temperatures (LnSR). Effect sizes above and below zero indicate a positive and negative response to warming, respectively. by a further 21% (p < .001) and 18% (p = .049), respectively. There was also a further 26% decline in performance of echinoderms from tropical (4-23°) latitudes, although this group still failed to exceed the threshold for significance (p > .050). Subtropical (26-35°) echinoderms were no longer significantly impacted by warming under the end-of-century scenario (p = .384); however, the percentage change in performance remained the same.

| DISCUSS ION
For echinoderms, finding patterns in vulnerability to warming be- Pereira Sampaio et al., 2021). A comparable response could have been anticipated for development success, as this would also be expected to shift due to an increase in the rate of cell processes (Munday et al., 2008;Pörtner et al., 2017;Pörtner & Peck, 2010;Schulte, 2015). However, there was no evidence that warming increased the speed at which early life stages reached developmental milestones, and more than ~3°C of warming was particularly detrimental to development success. Consequently, when considering warming consistent with end-of-century predictions, development success declined by 35%, indicating a possible negative effect of temperature on the cell cycle (Schulte, 2015), which may be further exacerbated by other anthropogenic stressors, such as ocean acidification Kroeker et al., 2010;Przeslawski et al., 2015). Developmental success may be more closely linked with survival, considering that slower development may indicate incompetency, and is associated with a higher risk of predation in the field (Munday et al., 2008;Ross et al., 2011).
Mortality due to supra-optimal thermal exposure will ultimately lead to declines in abundance, population viability, and even extinctions Eisenlord et al., 2016  . Surprisingly, although echinoderm survival declined with thermal challenge in the current analysis, limited impacts on other biological responses such as growth, feeding or movement were found. This is concerning, as previous studies have considered these metrics as indicators of sensitivity to climate change (Pereira Sampaio et al., 2021). Since these biological responses were not affected by warming in the models presented herein, echinoderms may maintain these behaviours and processes in a warmer ocean, leading to an energy deficit, and subsequent mortality. Behaviour is often expected to be the first response to environmental change as a means for individuals to buffer the negative effects (Wong & Candolin, 2015). However, since echinoderms are relatively limited in their mobility, they may have reduced capacity to alter their behaviour in response to warming. Given the relatively short experimental durations in most ex situ studies, it must be considered that behavioural responses to temperature may occur after lengthier exposure durations. Nevertheless, if echinoderm metabolic rate and hence energetic demands increase greatly with warming, while feeding and nutrition and therefore energy intake do not increase proportionally, an overall energy shortfall is expected during projected warming periods (Lawrence, 1984;Munday et al., 2008;Pörtner & Peck, 2010;Schulte, 2015). could not be explored here due to the lack of studies on brooding species (Menge, 1975;Pechenik, 1999 (Reich et al., 2015). The pattern may, however, be linked to variations in metabolic processes and the relative allocation of resulting energy to various body components, behaviours (e.g. movement) and processes (e.g. reproduction) within the body, which differs between taxonomic classes (Lawrence, 1984;Whitehill & Moran, 2012). A comparative study on echinoid (Arbacia punctulata) and ophiuroid (Ophiocoma alexandri) larvae, found that metabolic rates were lower in ophiuroids, suggested to be a consequence of their lower feeding rates which cannot sustain greater energy requirements (Whitehill & Moran, 2012). A lower metabolic rate in the ophiuroid larvae compared to the Echinoid larvae may be advantageous when food is limited, but a concomitant slower development rate may offset this benefit (Whitehill & Moran, 2012). In the present analysis, asteroids were deemed the most negatively affected by warming (30% decline in performance) and their performance declined by an even greater extent when only considering responses where warming magnitudes were in line with end-of-century predictions. In contrast to asteroid performance, ophiuroids showed trends of enhanced performance with warming.
However, due to the paucity of studies on this class, this analysis may not have had the resolution to detect any negative effects. An ability to make generalisations of thermal performance in different taxonomic groups may be helpful when estimating the responses of species, particularly those that have the ability to alter the structure and function of an ecosystem, and further research on underrepresented taxa are warranted (Ling, 2013;Menge et al., 2016).
The findings in this meta-analysis, to some extent, supports the general trend that the vulnerability of marine ectotherms to elevated temperatures increases with latitude from the poles to the equator, likely due to the increase in thermal stability (Hughes et al., 2018;Pinsky et al., 2019;Poloczanska et al., 2016;Sunday et al., 2011;Woolsey et al., 2015). When all data were included, echinoderms from subtropical latitudes, were the most negatively impacted by

| Data gaps, limitations and future directions
The present meta-analysis identified a number of data gaps and limi- These data gaps and limitations gave rise to the suggestion of several key considerations for future experimental studies on the impacts of ocean warming on echinoderms: • An increased effort should be made to record, not only coordinates of the collection location, but information regarding the habitat and the geographic range of the species.
• There should be a greater focus on studying groups of echinoderms that are poorly represented in the literature, as well as those that are of high ecological importance.
• There should be an accelerated effort in publishing studies that span multiple life stages and generations.
• All relevant data required for meta-analyses should be incorporated into the papers, including temperatures, means, standard deviations and sample sizes.
• All results should be published, even if they show non-significant effects.
A greater number of studies overall would provide the resolution to detect patterns in thermal responses with greater reliability.
Moreover, this would allow for the assessment of interactions between multiple predictors in meta-analytic models, to better forecast the vulnerability of echinoderm species into the future. To provide a more realistic assessment however, it is necessary to consider the myriad of other environmental and ecological factors (e.g. ocean acidification, food availability and predation) that may exacerbate, or even alleviate, echinoderm vulnerability to a warming ocean (Kroeker et al., 2010;Lucey et al., 2020;Melzner et al., 2020;Sampaio et al., 2021;. An important next step would be to incorporate multiple stressors associated with climate change, in a meta-analysis on echinoderm performance, that may have synergistic, additive or antagonistic effects when combined with warming (Harvey et al., 2013;Pandori & Sorte, 2018;Przeslawski et al., 2015;Sampaio et al., 2021).

| CON CLUS IONS
This meta-analysis provides a critical first step in understanding patterns of thermal responses among echinoderms and adds to the evergrowing body of literature predicting the negative impacts of climate change on marine species Dahlke et al., 2020;Pereira Santos et al., 2021;Poloczanska et al., 2016;Sampaio et al., 2021;Sunday et al., 2011). This taxon is clearly vulnerable to ocean warming, and many of the patterns observed in this meta-analysis fitted with broad expectations for marine ectotherms. For instance, proportionally greater negative responses to warming were observed in subtropical and tropical echinoderms and at larval ontogenetic life stages. However, there were some interesting patterns that did not fit initial expectations, including a lack of thermal vulnerability for many biological responses (i.e. growth, feeding and nutrition, movement, and reproduction).
The impacts of ocean warming and accompanying marine heatwaves are already evident in echinoderm populations, with diseaseinduced die-offs in numerous species being linked to elevated temperatures (Clemente et al., 2014;Dungan et al., 1982;Lester et al., 2007;Menge et al., 2016). Despite this growing threat, echinoderms may be able to persist in warmer oceans through range shifts and/or acclimation and adaptation (Ling et al., 2009;Munday et al., 2008;Poloczanska et al., 2016). However, there is a growing consensus that future ocean warming may be too rapid and marine heatwaves may be too abrupt for these mitigation strategies to keep

ACK N O WLE D G E M ENTS
The authors acknowledge the Traditional Owners of the land and sea country where the work was carried out, and pay respect to their Elders, both past and present. This work was supported by funding from the Australian Research Council (CE140100020;

FT190100015) and a James Cook University Postgraduate Research
Scholarship provided to BJL. The authors would also like to acknowledge the support provided to BJL by the Australian Research Council Centre of Excellence for Coral Reef Studies and the AIMS@JCU strategic partnership. Open access publishing facilitated by James Cook University, as part of the Wiley -James Cook University agreement via the Council of Australian University Librarians.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no competing interests.

A PPEN D I X 2
Data sources included in the meta-analysis. * indicates papers measuring metabolic rate that were removed from the main analyses for all moderators except 'biological response'. ** indicates papers where the experimental temperature was below the near-future mean annual temperature projection, and were therefore removed from the associated subset of data. A PPEN D I X 3 Predictors of variation in echinoderm vulnerability to warming, tested using separate multivariate linear mixed effects models. The predictor name and groups within them are provided.

A PPEN D I X 6
Sensitivity, publication bias and non-linearity tests methods and results.

Sensitivity analyses Methods
The presence of influential observations was tested by ranking the data by the effect size and removing the observations with the largest 10 effect sizes (in either direction) in a stepwise fashion and re-running all categorical models (i.e., biological response, ontogenetic life stage, taxonomic class, and region), to see if the overall significance of Q M , and the significance of individual groups of a predictor were affected.
The effect of pseudoreplication was then explored by removing publications one at a time that contributed two or more 'studies' due to measuring responses of multiple species and observing any changes in significance. Again, this was conducted for each of the categorical predictors. Lastly, the ontogenetic life stage and taxonomic class models were re-run, removing the three species that contributed the most data points: Acanthaster spp. (tropical asteroid; n = 110 data points), Apostichopus japonicus (temperate holothuroid; n = 78 data points) and Strongylocentrotus droebachiensis (temperate and polar echinoid; n = 53 data points).

Results
These models were robust to the removal of the 10 responses with the highest effect sizes (in either direction). Removal of papers that contribute two or more studies, in most cases did not change the model outputs to a significant extent. However, the responses reported in  which contributed six studies to the analysis appears to be particularly influential in two models. In the ontogenetic life stage A PPEN D I X 5 (Continued) model, gametes (and fertilisation) were no longer being significantly impacted by warming when responses from this paper were removed (p = .661). Furthermore, in the region model, echinoderms from tropical (4-22°) latitudes became significantly negatively affected by warming (p = .025), when  was removed. Moreover, removing  that contributed four studies to the meta-analysis, resulted in sub-tropical (26-35°) echinoderms being no longer significantly impacted by warming (p = 0.079), along with the overall Robust Q M (p = .077). Removal of Acanthaster spp. resulted in tropical echinoderms becoming significantly affected by warming (p = .040), and removal of Apostichopus japonicus resulted in the overall effect of the model becoming non-significant (p = .086). In all other instances the significances and direction of the response remained the same for the overall predictor models, and the individual groups.
In the tables, the groups in which the significance changed from the original models are in bold.
Life stage model output after removal of . Region model output after removal of .